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We consider Monte Carlo algorithms for the simulation of charged lattice gases with purely local 
dynamics. We study the mobility of particles as a function of temperature and show that the poor 
mobility of particles at low temperatures is due to "trails" or "strings" left behind after particle 
motion. We introduce modified updates which substantially improve the efficiency of the algorithm 
in this regime. 



I. INTRODUCTION 

The properties of many condensed matter systems can 
not be understood without considering the Coulombic in- 
teraction. DNA, proteins, polyelectrolytes, colloids and 
even water are all structured by electrostatics, which 
must be reproduced faithfully in any numerical study. 
Unfortunately, the simulation of the electrostatic inter- 
actions is difficult; due to the slow decay of the potential 
in 1 jr one can not truncate the interaction ^ as is often 
done with other molecular interactions. Most working 
codes now use a variant of the Ewald sum to account for 
the interaction between periodic images of the basic sim- 
ulation cell. As a consequence the time need to evaluate 
the electrostatic interaction can dominate in the simula- 
tion of charged systems. 

In a Monte Carlo simulation with the Ewald method, 
the motion of a single charge requires summing its in- 
teractions with the TV — 1 other charges and their pe- 
riodic images, resulting in a O (TV 2 ) computational cost 
per sweep. This becomes impractical when TV is large. In 
simulations with explicit modeling of all charges TV > 10 4 
is commonly required. In molecular dynamics the situ- 
ation is better: when all particles are moved simultane- 
ously, better CPU time scalings are possible ranging from 
0(N) for multigrid algorithms to (9 (TV 3 / 2 ) for an op- 
timized Ewald summation 0. However these molecular 
dynamics codes are complex to implement. 

The unfavorable complexity of conventional Monte 
Carlo methods originates in the use of the electrostatic 
potential which is the solution to Poisson's equation 



-p/e. 



This equation has a unique solution for given charge dis- 
tribution and boundary conditions. When a charge is 
moved, the new solution for <E> is computed and the inter- 
action energy g^$(r^) of the moved charge with all other 
charges qi in the system changes. The electrostatic inter- 
action implemented in this way is instantaneous. Note 
however [j, the thermodynamical study of charged sys- 
tems does not require instantaneous Coulombic interac- 
tions: The free energy is also correctly sampled when 
only Gauss's law 

divE p/e 



is imposed on the electric field. The fact that solutions 
to Gauss's law are not unique results in an extra flexibil- 
ity which allows one to implement a purely local Monte 
Carlo scheme for the simulation of systems with electro- 
static interactions. The computation effort is reduced 
to 0(N) per Monte Carlo sweep. The disadvantage of 
the algorithm is that it requires a grid to discretize the 
electrostatic degrees of freedom, however this is also true 
of multigrid and Fourier methods used for molecular dy- 
namics. 

The final efficiency of the Monte Carlo algorithm de- 
pends on the number of sweeps required to sample inde- 
pendent configurations which in turn is a function of the 
particle mobility resulting from the Monte Carlo dynam- 
ics. Highly mobile charges enable one to generate inde- 
pendent configurations rapidly; if charges were to become 
"trapped" or "localized" due to their interaction with the 
field it could prevent the generation of uncorrelated sam- 
ples. Monitoring the acceptance rate of particle updates 
may only give partial information in that on the total 
efficiency of an algorithm. For instance trapped particles 
could move locally (resulting in a good acceptance rate) 
without being able to explore all of space. 

In this article we perform a detailed study of the charge 
mobility p in local Monte Carlo algorithms in order to 
compare efficiencies of various implementations. Firstly 
we develop a technique to measure p by relating the mo- 
bility to the dynamics of the average electric field, E. 
The mobility will be studied as a function of tempera- 
ture. With our previous implementation, p drops dra- 
matically at low temperatures, becoming unmeasurable 
for parameters which are needed to study typical ma- 
terials: For instance monovalent ions in water at room 
temperature where e = 78, T = 300 K, a = 1 A with 
a the mesh size. The drop in efficiency originates in the 
constrained dynamics of the electric field, leading to the 
generation of "trails" or "strings" which trap particles at 
low temperature and suppress their mobility. 

We will explore ways of reducing this trapping. The 
update law introduced previously 4( for particle motion 
is not the only way one can move a charge. Even if 
each charge update must be accompanied by some field 
update, the latter is only loosely constrained. Duncan, 
Sedgewick, and Coalson recently used this fact to in- 
troduce p a better particle-plaquette update. We will 
present several field updating schemes leading to less 
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trapping. These schemes are very flexible in that they 
have a freely adjustable "spreading" parameter w, upon 
which their effects and their computational complexity 
depend. Schemes which use a larger spreading param- 
eter are more time-consuming but lead to much larger 
efficiencies at physically interesting temperatures. 

We have already shown 6] that off-lattice implemen- 
tations of the algorithm (using continuous interpolation 
of charges with splines) do not suffer the mobility drop 
that we discuss in this paper. Rather, our present work 
is motivated by the existence of a wide spectrum of inter- 
esting and useful lattice models. For example, it is known 
0, 0, E3 that finely discretized lattice fluids exhibit the 
same critical behavior as continuum fluids. Another ex- 
ample is the bond fluctuation model for polymers [TH HH 
which is rather easily generalized to study charged poly- 
mers, or polyelectrolytes. All these models already use 
"spread" or extended particles where the hard cores of 
the particles span several lattice sites in order to reduce 
lattice artefacts to an acceptable level. 

We will begin (Sec. [DJ with a description of the the- 
oretical basis and implementation of local Monte Carlo 
algorithms with electrostatics. In Sec. IIIII we show how 
to measure the mobility of charges and apply the method 
to the simplest algorithm. We interpret the behavior of 
the acceptance rate and mobility as a function of tem- 
perature fSec. HV|) . introducing the concept of field trails 
or strings. We show how to increase particle mobility in 
Secs. lyllVll Finally we will present the CPU time for rep- 
resentative simulations and give the reader an estimate 
of optimal parameters. 



II. MONTE CARLO ALGORITHM 



We give a basic description of the algorithm previously 
developed in Refs. @,[l3,u3- We first recall its theoreti- 
cal basis, and then present the simplest implementation, 
highlighting places where the method can be further op- 
timized for speed. 



A. Theoretical foundation 

In order to sample configurations of a system contain- 
ing charges, we require that a set {r^} of particle positions 
is generated with weight 

z({r z }) = Jy^WWr;{r.}) d 3 r + u({^})} 

where p(r; {r^}) = qid{r — r^) is the charge distribu- 
tion of the configuration, <l>(r; {r^}) is the unique solution 
to Poisson's equation, and ^/({r^}) is the potential of all 
other interactions. The partition function then is 



In the usual treatment of electrostatic interactions the 
electric field is given by Ep = — grad^>: it is unique and 
satisfies both Gauss's law divEp = p/e and the static 
version of Faraday's law curlEp = 0. We chose the 
convention where the potential of a charge q is q/Airer. 
The algorithm is based on relaxing Faraday's law so that 
E = Ep + curlQ, a decomposition familiar from the 
Coulomb gauge of electrodynamics. Fourier transform- 
ing we find that Ep is longitudinal: 

(Ep) k = -ik$ k || k 

and that curl Q is transverse: 

(curlQ)k = ikxQ k _L k. 

As a consequence, the electrostatic energy 

e - J E 2 d 3 r=^(y Ep 2 d 3 r + J (curlQ) 2 d 3 r) 
= \ [ p$d 3 r+ J f (curlQ) 2 d 3 r 

2 JV ^ Jy 

so that the statistical weight of a configuration of charges 
and field is 

z'({r l },E ± )=e-^ U{{r ' })+1 2lv E2d!ir } 



; ({r0)e" /3 f Iv 



where for clarity we have introduced the transverse field 
E^ = curlQ. This field is constrained by divE^ = 0; 
it is independent of the charge configuration. Thus, the 
partition function of the system of charges and field splits 
into two parts: 

Z'= /(J] d 3 r^VE ± 8(divE ± )z'({r l },E ± ) 

- (/(n <fr<)««'<») a) 

= Z X >Et r 

where Z tY is the partition function of transverse field. 
The statistical weight of a configuration of charges is 
z({r^}) x Z tr ; all the weights have been multiplied by 
the same constant. Hence configurational probabilities 
are left unchanged. Of course sampling this system re- 
quires introducing Monte Carlo moves appropriate for 
integrating over E^ degrees of freedom. 



B. A charged lattice gas 

We consider a cubic simulation cell of L 3 sites with pe- 
riodic boundary conditions. Particles are placed on sites 
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of a lattice with mesh spacing a, and field variables rep- 
resenting electric flux are defined on links. The electric 
field divergence at a site is the sum of fluxes over the six 
outgoing links. (See Fig. QJ) Transverse field degrees of 
freedom appear as a nonzero line integral of E on pla- 
quettes of the lattice. 





FIG. 1: Left: cubic lattice mesh around site i. The elec- 
tric flux faj flowing upward through the hashed cube face 
is assigned to link on the right: Eij — faj/a 2 . Field 

divergence at i equals the outward flux through the cube sur- 
face, X^g{nn} Ei,j where NN stands for nearest neighbors. 
In following figures only two dimensions of the lattice will be 
shown for clarity. 

To start a simulation we must construct a state con- 
sistent with Gauss's law. We initialize the electric field 
for the simulation with a single sweep through the net- 
work: We use a procedure that follows a Hamiltonian 
path through the lattice. Such a path visits each site just 
once and traverses each link either once or zero times. 
We begin by initializing all field values on the lattice to 
zero and start at an arbitrary point 1 of the lattice; the 
node 1 holds the charge q±. A single link of the path, 
{1,2}, connects it to site 2, on which we set the outgoing 
field to qi/a 2 e; Gauss's law is now fulfilled on site 1 and 
we move to the node 2. 

At each step, on arriving at site i holding q^ we have 
already solved the Gauss constraint for sites {1, . . . , i— 1}. 
The incoming link to the site, {i — l,i}, thus bears 

the initialized field Ei-ij = Qj^ / a<2e - We 

now set the outgoing field to ^X}j'=i /& 2 e so 

that Eij-i + Eij+i = qi/a 2 e: Gauss's law is now ful- 
filled on site i and we go to site i + 1. At the end 
of the path, we reach site V = L 3 with Ey-iy = 

(X^=:T _1 /& 2 e. The imposition of periodic bound- 
ary conditions in charged systems is only possible if the 
total charge Q is zero (otherwise the total energy is di- 
vergent). Thus Ey-iy = (Q — qv)/a 2 £ = —qv/^e and 
Gauss's law is satisfied everywhere on the lattice. 

We take advantage of the new field degrees of freedom 
to construct local updates for charge moves. Consider 
an initial configuration p^ where a charge q is at point 
A, and the initial electric field E^ satisfies Gauss's law 
(divE^ = pW/e). If a trial places q at point £>, the final 



configuration is p^ = p^ — g£(r — r^) + q5(r — re), and 
a new solution for the field must be found. In order to 
remain consistent with Gauss's law it is sufficient to add 
to field lines £E flowing from B to A and totalizing 
q/e flux. The final field E^ f ) = E« + 5E now satisfies 
again Gauss's law for the final charge configuration: 



,(0 



d ivE« = 9 — + U{t - r B ) ~ -S(r - r A ) 



We call 5E "the slaved update". Nothing having been 
required of it except the total flux, we can choose it to 
be localized in space, so that charge moves result in local 
updates of the simulated system. In addition 5~E should 
be symmetrically chosen so that detailed balance applies 
to the forward and reverse updates. Our previous choice 
of (5E, which modifies just the link connecting A to B is 
illustrated in Fig. Nevertheless, the "total flux" con- 
straint lets us free to use more complex slaved updates. 
We will show in this paper that splitting £E into several 
lines helps increase the algorithm efficiency. 
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FIG. 2: A pair of lattice sites, before and after a particle 
move. Left, the initial configuration is made up of a charge 
at A and one solution to Gauss's constraint: at A the field 
divergence is six times 1/6, equaling q = 1 (in reduced units 
where a = 1, e = 1), and at B it is 5 x 1/30 - 1/6 = 0. Right, 
the charge has moved to B and a flux SE = q = 1 flowing 
from B to A has been added to the central link. Then Gauss's 
law is again verified: at A, divE = 5x1/6 — 5/6 = 0, and at 
B, divE = 5 x 1/30 + 5/6 = 1 = q. 

Finally in order to correctly sample the partition func- 
tion we integrate over the transverse degrees of free- 
dom of the electric field. We do this with Monte Carlo 
moves which change the circulation of the electric field, 
but do not modify its divergence. One way of doing this 
is by modifying the field on the four links defining a pla- 
quette [Fig.|JJfa)]. If one increases the field on links along 
the edge of a given plaquette by some constant value, at 
all sites divE remains constant. This kind of update, 
being local, leads to diffusive dynamics for E^: 0(L 2 ) 
sweeps are needed to yield an independent configuration. 

An alternative method of integrating over E^ was in- 
troduced 14]. These worm updates [Fig.Efb)] make use 
of a biased random walk to generate a closed contour 
along which the field is modified. This contour visits typ- 
ically L 3 sites and turns out to be particularly efficient 
at equilibrating the electric field at all length scales si- 
multaneously: all Fourier modes of E^ decay at the same 
rate, 0(l) sweeps are enough to produce an independent 
field configuration. 
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FIG. 3: To update one may choose between plaquette 
moves (a) which increase the field along a single plaquette 
edge, and worm moves [(b), dashed line] which modify it along 
the path of a given random walk. Both are thermodynami- 
cally equivalent since any field configuration reached through 
worm moves may be obtained through multiple plaquette up- 
dates [as shown in (b)] combined with updates of the k — 
mode of the field. 



The aim of a Monte Carlo algorithm is to produce sta- 
tistically independent configurations with minimum com- 
putational cost. The local updates described above allow 
one to efficiently update charge and field configurations. 
However in order to understand the global dynamics and 
convergence of the algorithm we shall study electric field 
autocorrelation functions. We now show that high mo- 
bility p of the charges leads to fast decay of the field 
correlations. 



III. MEASURING CHARGE MOBILITY 

Under the dynamics of the algorithm, Figs. and E 
remains consistent with Gauss's law at all times. Con- 
sidering the time derivative of this law, we find 



div 



dE <9(divE) 



I dp 

lot 



dt dt 
which translates in the implementation as 

1 



div 5E ■ 



-Sp. 



(2) 



Updates to the electric field can be considered as being 
due to local currents such that 



div J 



dp 
St 



o, 



where we introduced the time unit St = 
step. Combining Eqs. (0) and (0) we find 

divJE = — — div J, 

e 



(3) 

1 Monte Carlo 



or 



SE 



St 



J + 5t curlH 



(4) 



with H arbitrary. This is a discrete version of Ampere's 
law of electromagnetism. 

Spatially averaging Eq. (gj), we find the change in the 
average electric field E during an update, 



SE 



St 



J. 



(5) 



This equation is independent of H; the last term in 
Eq. (gj) gives zero due to periodic boundary conditions. 
This is consistent with the fact that local plaquette up- 
dates do not change the average electric field in a periodic 
system. 

Our simulations are on a system containing N mo- 
bile unit charges, either the symmetric plasma made up 
with N/2 particles of each sign = ±e), or the one- 
component plasma (OCP) of N positive charges moving 
in a fixed negative background. Linear response gives in- 
sight on the relation between charge mobility and field 
evolution. The electric current is due to the movement 
of mobile charges, 



(6) 



where i £ {+,—}, pi are charge densities and rii are num- 
ber densities; n_ = for the OCP. On average, velocities 
are related to field by 



pq%E. 



(7) 



Given the charge symmetry of the algorithm, positive 
and negative ions in a symmetric plasma have the same 
mobility. Equations (0 and J7J) lead to J = e 2 (n + + 
n-)pE. n + + n_ = n 
mobile charges. Hence 



N/V is the number density of 



J e 2 npE. 



(8) 



We should bear in mind that these relations are phe- 
nomenological. For example, in Eq. Q proportionality 
holds only when the field intensity is not too high. It will 
also become apparent that in certain limits p can fall to 
zero for large, dilute systems. 

Substituting Eq. JSJ) in Eq. JJj, and replacing the dif- 
ference equation by a differential equation we find that 



dE 

~dt 



e 2 np 



E. 



(9) 



E is the k = Fourier mode of the electric field. Equa- 
tion ® implies that the autocorrelation function of this 
mode behaves as follows: 

_ _ e 2 n/i 

(E(t')E(t' + t)) t ,=Ce-— *, 

where C is the squared amplitude of the thermal fluc- 
tuations of E. Measuring this autocorrelation function 
we find exponential decay with a characteristic time 
tq = e/(e 2 np) or equivalently a decay rate Aq = e 2 np/e. 
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We fit all our numerical data with a single exponential 
and verify the quality of the resulting curve by eye. 

In our simulations we also monitored other modes of 
the field and found that the mode k = is the slowest. 
Higher modes of the field couple directly to plaquette 
updates as well as particle motion, and relax with the 
dispersion law = Ao + Dsk 2 [l2|. Larger k are less 
sensitive to low particle mobility (low Ao). They will not 
be considered further in this paper. 

The time scale To can be understood with a scaling 
argument. In order to produce two uncorrelated sam- 
ples of the system, one should wait for the charges to 
diffuse through the characteristic correlation length of 
the system, the Debye length Id — \feksT / e 2 n. Thus 
to = Id 2 /D, and Ao = kd 2 D = e 2 nD / 'eksT ', with a 
diffusion constant D. We recover the above expression 
for the relaxation rate if we use the relation D = ksTfi. 
D defined in this way relates to the mobility of charges 
under an external electric field (where opposite charges 
move opposite ways), not to the mobility under a con- 
stant external force (where all particles move together). 

Thus we will measure the mobility of charges or, equiv- 
alently, their diffusion coefficient, by computing the au- 
tocorrelation function of the average electric field. This 
method has the additional advantage that we need not 
keep track of the winding of particles across the periodic 
cell boundaries. 



IV. LIMITING FACTORS FOR MOBILITY 



Acceptance rate is often used to monitor the efficiency 
of Monte Carlo simulations. However a high acceptance 
rate does not necessarily mean useful work has been per- 
formed. For example, the diffusion dynamics of point 
defects or interstitials in a crystal are very slow. In 
a Monte Carlo simulation most time is spent vibrating 
the atoms around their equilibrium positions; even in 
the limit of very rare diffusion events the acceptance 
rate of trial moves remains appreciable. Another ex- 
ample is magnetization reversal of an Ising ferromagnet. 
The state where all spins are oriented against an applied 
field is metastable but with very long lifetime; since the 
Metropolis algorithm is already very inefficient one uses 
rejection-free algorithms [13 to update individual spins, 
unfortunately after a spin has been flipped it is almost 
certainly flipped back at the next step, so that the magne- 
tization never reverses within accessible simulation times. 

With our algorithm, a simulation performed at very 
low density, Fig. EI shows that the diffusion coefficient D 
of charges drops much faster at low temperatures than 
the acceptance rate of particle moves: atoms simply wan- 
der around their mean locations, rather like in the exam- 
ples above. 
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FIG. 4: Logarithm of acceptance rate R, and of diffusion 
coefficient D expressed in a 2 per particle sweep, versus inverse 
temperature. Solid line, Eq.[H3 dotted line, guide to the eye. 
R and D are close to 1 at high temperatures, but D drops 
much faster than R on decreasing T. One component plasma 
of two positive unit charges, box of size L — 15. 



A. Variation of the acceptance rate 

In the algorithm summarized in Sec. Ill Bl motion of a 
charge modifies the field on the single link along which 
the particle has moved (see Fig. 0). Let Eab be the 
field intensity on this link before the move. The diver- 
gence at A is g/a 2 e, and since in the absence of other 
nearby charges E must be isotropic around A we expect 
that {Eab) = q/6a 2 e. Fluctuations of imply that 
Eab — q/Qa 2 e + 77, with 77 a Gaussian random variable 



with standard deviation a. The energy in these fluctu- 
ations is 3Z/ 3 a 3 (f77 2 ) = §eL 3 a 3 a 2 . From equipartition 
and given that there are two polarizations of E_l, the 
energy in Ex is also approximately L 3 ksT, thus we con- 
clude that a 2 = 2k B T/3a 3 e. 

During motion of the charge, the field on AB is modi- 
fied to — 5q/6a 2 e + 77. The energy difference between the 
two configurations is thus SS = qa(q/3a 2 e — 77). With 
the Metropolis algorithm when 77 < q/3a 2 e the trial is 
accepted with probability exp(— SS /fe^T), otherwise it is 
automatically accepted. Computing the average over all 
values of 77, we find the acceptance rate 



R = erfc (q/2^3aek B T s j . 



(10) 



We plot this function together with numerical results in 
Fig. When T is small, the asymptotic expansion of 
erfc gives 



R 



12T 



e-^ 12T (l + 0(T)). 



(11) 



Defining y = In (R/y/T} and x = 1/T, we find an Arrhe- 
nius law for the acceptance rate y = —x/12 + const. + 
G(l/x), which is illustrated in the inset of Fig.0 
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FIG. 5: Acceptance rate R of charge moves versus tem- 
perature. + , simulation results; solid line, Eq. EH Inset, 
y = \n(R/VT) against x = 1/T. Numerical results approach 
the asymptote Eq. (fTTj) of slope —1/12 (dashed line). A pair 
of opposite charges, box size L — 15. 



We conclude that particle motion becomes hard with 
this update scheme for T < 1/12 due to a finite energy 
barrier. One of our aims in the rest of this paper will be 
to reduce the barrier so that the acceptance rate remains 
high even for temperatures T <C 1/12. 



B. Field trails and string tension 

Let us consider two closely separated charges with the 
electric field in equilibrium. The field has the usual dipo- 
lar form familiar from elementary electrostatics. What 
happens if we pull very hard on the positive charge so 
that the separation between the charges increases rapidly, 
without updating the plaquette degrees of freedom? Dur- 
ing the motion of the particle each link traversed is mod- 
ified by q/a 2 e leaving behind a "trail" of modified links 
(Fig.EJ). With time the field configuration will relax back 
to a dipolar form because of the updates of the plaque- 




FIG. 6: The field update produced by successive moves of a 
charge q is a field string of intensity SE — q/a 2 e connecting 
the particle to its starting position (dashed circle). The trail 
remains as long as no plaquette update intervenes. 



ttes equilibrating the transverse field. However on a short 
time scale there are few plaquette updates, and dragging 
the charge along r links costs an energy which we can 
estimate to be r7o, where 



70 q 2 /2ae 
is our estimate of the energy rise per link, 



(12) 



-aqE 



2ae 



and E has zero mean. There is a "string tension", 70, 
pulling the particles back. 

In the presence of an external electric field, a pair of 
opposite charges normally separates. A finite string ten- 
sion implies that this mobility is suppressed. One must 
spend much numerical effort on updating the plaquettes 
in order to destroy the trail and stop particles backtrack- 
ing. While the string tension is positive at low tempera- 
tures we will now argue that the thermodynamic tension 
7 should become zero at a finite temperature: above it 
the mobility is high even at low frequencies of updates in 
the plaquettes. 

Consider a trail joining two fixed test charges separated 
by a distance r and let the length of the trail joining 
them be £. If £ ^> r we can estimate the number of such 
paths from the statistics of the path: Nz = 0(z^. z 
is a connectivity constant characterizing the geometry of 
the walk. For a random walk z = 6, for a self- avoiding 
walk z = 4.68. We now estimate the free energy of the 
configuration as 



F^ijo- k B T£lnz = j£. 



(13) 



From this expression we can expect two distinct dynamic 
regimes for the algorithm. At low temperatures the ten- 
sion 7 is positive and it is most favorable for the trail to 
remain short, £ ~ r. The free energy for separating the 
charges is indeed linear and we have a phenomenon simi- 
lar to confinement in gauge theories. This confinement is 
only destroyed by the dynamics of the plaquettes which 
slowly relaxes the trail into a dipolar field configuration. 
At temperatures higher than T c ~ jo/ks^z ~ 0.3 the 
line tension drops to zero and the particles become un- 
confined. Even without plaquette updates the particles 
remain mobile and can separate easily. 

We also note that there is a very close analogy be- 
tween this picture of roughening trails and the (2 + 1)- 
dimensional Hubbard model in the phase approximation, 
which can be expressed as a set of fluxes on a lattice 
0,^3- This model has two thermodynamic phases, one 
with tense field lines which are strongly suppressed, and 
a superconducting phase in which field lines proliferate. 
The transition occurs at a temperature T « 0.33. 

In Fig. 21 we used a split in which one half of all up- 
dates try to move one of the two particles, and one half of 
updates modify a randomly chosen plaquette. The num- 
ber of plaquettes (3L 3 « 10 4 ) is much larger than the 



7 



number of particles (TV = 2), so that a given plaquette 
is rarely updated, trail formation is probable. The dif- 
fusion coefficient of charges indeed drops at a crossover 
temperature T c w 0.2 which qualitatively agrees with the 
above estimate. 

How do we expect this trail-limited mobility to vary 
as a function of charge density? If a charge i creates a 
trail, and a charge j of the same sign crosses it, then 
j will also feel the mean force mentioned above. If j 
is now dragged back along the track of z, the field up- 
dates will erase the trail (Fig. Q). Afterwards, neither i 
nor j are linked to their initial positions. We thus ex- 
pect that the effect of the trails is cut off at a distance 
comparable to the inter-particle spacing. Indeed we do 
find that the mobility increases on simulating systems of 
increasing charge densities. Thus in this paper we will 
concentrate on improving the efficiency of the algorithm 
at very low densities, working most often with samples 
containing just two charges. 




FIG. 7: A second charge joins the string left by the charge 
of Fig. El it is dragged along the original path. Field updates 
then erase the previous trail (dashed line). The field strings 
no longer connect the particles to their respective starting 
sites (dashed circles). 

In the next two sections we will modify the slaved 
updates in order to reduce the bare tension of strings. 
With a lower 70 we will lessen the crossover temperature 
T c w 0.2 which results from the balance between energy 
and entropy expressed in Eq. ([T3|) . To efficiently simulate 
condensed matter systems, particles must remain highly 
mobile to much lower temperatures: T = 0.2 corresponds 
to Ib = a/4irT & 0.4 A (with a = 1 A), whereas the Bjer- 
rum length in water at room temperature is Ib ~ 7 A. 
Therefore we aim at lowering T c by a factor of approxi- 
mately 20. 



V. EXTENDED CHARGES 

The expression O for the bare tension of the string 
is quadratic in the charge, 70 = q 2 /2ae. Let us now split 
the string between two particles into K substrings; each 
substring carries a flux of q/Ke. The bare tension of 
each substring will be 70 / K 2 , and that of the whole split 
string will be K(jo/K 2 ) = jo/K. 



In this section, to form split strings we spread the par- 
ticles on cubes of side w; each site in the cube carries a 
subcharge of q/w 3 , and when a particle moves the field is 
updated on the w 3 links crossed by each subcharge. We 
use values of w ranging from 1 (the original algorithm) to 
5, and measure the acceptance rate of particle updates. 
When we plot the rate as a function of w 3 T (Fig. we 
find that all curves collapse, except at low temperatures 
for the two opposite charges due to pairing. We also sim- 
ulated point charges with the coupled update proposed 
by Duncan, Sedgewick, and Coalson in Ref. Q (hereafter 
denoted by "DSC"), and the acceptance rates collapse 
equally well. 
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FIG. 8: Acceptance rate of particle moves versus temperature 
for OCP (open symbols, w = 1 to 5) and a pair of charges 
(filled symbols, w — 1 to 3), temperature rescaled by w 3 
(DSC, 7). 

The scaling of the acceptance rate in w 3 can be un- 
derstood as follows: motion of each part of the particle 
is hindered by a barrier which varies as (1/w 3 ) 2 . The 
barriers are additive leading to a local barrier with an 
amplitude which varies as 1/w 3 . 

When we plot mobility (determined from the dynamics 
of E) as a function of temperature, we find that the ben- 
efit obtained from charge spreading is not proportional 
to w 3 ; curves collapse on using a scaling with w 2 , Fig. [SJ 
The cross-section area of the extended charges is equal 
to w 2 , so their field trails are made up from K = w 2 field 
lines of strength q/a 2 w 2 e. This gives a bare tension for 
the trail of jo/w 2 . When trail formation limits mobility, 
the typical crossover temperature T c thus scales as 1/w 2 . 
This seems to indicate that the statistics of the paths and 
the connectivity constant do not change with w. 

To further confirm the idea that field trails are limiting 
mobility we introduced a new kind of field update: We 
define a cubic box of side b centered on a site occupied by 
a particle, and then generate a worm update (Sec. Ill Bl 
and Ref. [l4|) inscribed in the box. At each Monte Carlo 
step, the algorithm attempts one of three updates, either 
a particle move, a plaquette update, or a "local worm" . 
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FIG. 9: Mobility in OCP versus temperature (N = 2, 
L = 15). Each charge is spread on a ^-site-side cube. Data 
collapse when temperature and Ao are scaled by w 2 . For 
T — > oo the diffusion coefficient of particles is bound by 
1 a 2 (sweep) -1 : D saturates and Ao ~ 1/T. 



FIG. 10: Mobility of OCP particles versus temperature at 
low density (N = 2, L = 15). Here local worm moves are 
introduced (see text), and now data are made collapse when 
scaled by w 3 , showing that field trails have been removed, as 
opposed to data of Fig. 



Since all choices are reversible detailed balance is verified. 

Worm moves are known to lead to fast relaxation of the 
field, so these new "local worm updates" should allow one 
to spread out the field trails efficiently, concentrating the 
computational effort around charges, where the trails are 
formed. By introducing them in a 1:1 proportion with 
particle moves (with b satisfying b 2 > w 3 ), we expect 
to cancel the effective string tension. This computation 
is very expensive; one "local worm update" is far more 
costly than one particle update. We did not seek fur- 
ther optimization, and do not recommend this method 
for production of data with the algorithm. 

We find that the crossover temperature T c of the mo- 
bility drop decreases with these local worm moves. In 
Fig. the data superimpose if we rescale by a factor 
of w 3 , implying that the trails are no longer dominating 
the dynamics. The w 3 scaling may be indicating that 
dynamics are now limited by the local barrier to particle 
hops described in Sec. II V Al 

The spreading of particles over several sites clearly 
modifies the interactions at short distance. One should 
introduce a hard core interaction for distances less than 
wa, corresponding to the diameter ro of the particles. 
Much of the interesting physics in soft condensed mat- 
ter depends on the ratio of the Bjerrum length lg = 
e 2 /47re/c#T to the particle size. One is typically inter- 
ested in the range 5 < Z#/ro < 20, which corresponds 
to 0.004 < wT < 0.02. While we have succeeded in re- 
ducing the crossover temperature T c by a factor 1/w 2 we 
have also changed the physical length scale by a factor 
w. The final result is only a factor w improvement in 
T c when measured in physical units; a lattice algorithm 
suitable for condensed matter simulation would require 
w ~ 20. Such fine discretization has been used in lattice 
models to reproduce correctly thermodynamical proper- 



ties of some systems [l^l • However for cases where 
this is not required, one might prefer to avoid such large 
w. We now explore methods of moving charges which 
do not require permanent spreading so that the effective 
length scale in the simulation is not modified. 



VI. TEMPORARY CHARGE SPREADING 

There is a direct way of reconciling the requirements 
that charges are extended during their motion but oth- 
erwise point like: Before moving a particle, one should 
first spread its charge evenly onto neighboring sites, then 
move all sub chare; es as a block, and finally bring them 
back together (see Fig. [HJ- This defines a charge move 
involving three substeps. 

Each step consists of a set of currents. When a charge 
is split a current j^ 1 ^ flows from the central site. Motion of 
the particle generates a current, j( 2 ). When the charge is 
collapsed to a point a current j^ 3 ^ flows from the neighbor- 
ing sites back to the center. To maintain the constraint of 
Gauss's law, each of these currents is associated with 
a field update SE^ = -j^St/e. For step 2, the current 
on each modified link is = q/w 3 a 2 St, as above. Dur- 
ing step 1 the values of the current on links {i} are 
under-determined, they are constrained only by charge 
conservation Eq. (|3J). We thus additionally require that 

X^Of 1 ^) /2 be a minimum, giving a unique, reversible 
recipe for the current. We solve for j^ 1 ^ by minimizing 
the functional 
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FIG. 11: Temporary spreading of charge (here with w = 3). 
First step, spread the charge to w 3 sites. Second step, move 
the extended set of charges. Third step, inverse of the first; 
charge fractions collapse back to a point. The overall field 
update is the sum of the individual steps. 



The current is the solution of a Poisson-type problem, 

5pM 



-gradA and V A 



St 



with j^ 1 ) • n = on the boundary of the spread charge. 
The solution to this equation is computed once during 
initialization of the simulation and stored in a lookup 
table. Step 3 is the exact reverse of step 1, j( 3 ) = — 
On adding the fields jW, ]^ 2 \ and j( 3 ) we find a flow 
going from the starting site to the final site, and taking 
several paths. 

If now we simulate our test system using this version 
of the algorithm and plot the mobility of particles ver- 
sus temperature (Fig. IT2|) . we find practically the same 
results as in Fig. The crossover temperature scales 
with 1/w 2 . The advantage is that the particles are still 
pointlike unlike Sec.|3 Thus we have improved the low- 
est temperatures efficiently accessible by a factor of w 2 
without changing the physically important length scale. 

The method has the advantage of both simplicity and 
generality. 

• A small Poisson equation is solved once before 
starting simulation in order to determine the "cur- 
rent map". 

• One is free to choose as an intermediate state an 
arbitrary charge cloud. 

We note that temporary spreading includes the Monte 
Carlo algorithm of DSC as a special case, it is sufficient 
to consider the six nearest neighbors of a site, plus the 
site itself, as the volume over which the charge is to be 
spread. Each site thus gets one seventh of the total 
charge which explains the scaling of acceptance rate 
in Fig. El The result of steps 1 through 3 of Fig. El 



FIG. 12: Mobility of temporarily spread charges versus tem- 
perature. OCP, N = 2, L = 15. The mobility with no spread- 
ing (w = 1) is also given. Data are rescaled by w 2 (DSC, 
3.77). 



yields a current 3q/7a 2 St on the center link and q/7a 2 St 
around the four plaquettes adjacent to it. The curve ob- 
tained by using DSC updates collapses with the others in 
Fig. El when scaled by 49/13 = 3.77 which comes from 
a simple estimate of the bare string tension. However in 
our simulations we have not implemented the additional 
step of simulating the update with heat bath rather than 
Metropolis update. 

Rather than performing three successive steps, another 
implementation of the temporary spreading consists in 
precomputing the total "map" of currents j^ 1 ) +j^ 2 ^ 
Such an implementation is faster, avoiding multiple up- 
dates of the same links through steps 1-3. However, keep- 
ing the three sub-steps distinct allows one to perform 
several intermediate updates of step 2, moving the par- 
ticle a large distance before "recondensing" it to a point, 
as follows. Consider a (starting) configuration Cs of the 
system. We randomly choose one particle and spread its 
charge unconditionally; the field is updated accordingly 
and we label the new configuration as Ci=$\ the energy 
difference between Cs and Co is stored. Then we succes- 
sively try d moves of that particle in random directions, 
which lead to configurations C;, 1 ^ i ^ d. The field 
is updated and trials are accepted with the Metropolis 
probability 



m(A£) = min ( l,exp 



AS \ 



k B T 



)) 



where AS is the energy difference between the tried and 
current configurations. Finally the charge is condensed, 
yielding the ending state Ce, and the whole update is 
accepted with probability 



Pa 



i(AS(Cs^C ) + AS(C d ^C E )). 



The method saves computational effort, for each series of 
d moves there are only one spreading and one condensa- 
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tion steps; there would be d such steps if the procedure 
of local hopping of Fig. ITT1 were used. 

To prove that detailed balance is obtained, consider an 
instance of such an update. Its global probability is 

d-l 

P(C S - C E ) — Pace 

i=0 

The probability of the (i + l)th step is 
p(C l ^C l+1 ) = 

— m(A£(C^ — > Ci+i)) if Ci 7^ Ci+i (accepted trial), 
1 ~ \ Yl m ( A£ ( C i ^)) if C i = C i+i (rejected trial), 

ol' 

where the sum runs over the six directions of space, and 

A£(Ci —>) is the energy change corresponding to a trial 
move in the a' direction from configuration i. The prob- 
ability of the reverse update reads 

d-l 

Kc E ^Cs)=p / acc n^+i^c,), 

i=0 

where the global acceptance probability is 

p' acc = m(A£(C E - C d ) + Af (Co - C s )). 
When Ci = C i+1 , 

p(C,^C,+i) _ 1 f Af^+i' 



when ^ 7^ C*+i, 
p(C, C <+ i) 



and 

Pace 
Pace 

so that 

p(C s ^ C E ) 
p(C E - C s ) 



Ci) 



exp 



+1 



exp 



1 = exp 

m(Ag^ i+ i) 
" m(A5 <+ i^<) " ~^ V feT 

A£(C S -> Co) + Ag(C d -> C E ) 



A^s^o + E^^+i + Af d ^ E 



exp 



exp 



£(Ce)-£(C s ) 
k B T 



To check that the mobility is not changed when these 
long-ranged particle moves are used, we simulated our 
test system with them, fixing d = 15. Regarding time 
units, one such update amounts to d elementary Monte 
Carlo steps. On Fig. [E3 we find that the mobility of 
charges is very little affected by the use of long distance 
particle updates. However CPU cost is reduced, as is 
shown in next section. 
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FIG. 13: Mobility of temporarily spread charges versus tem- 
perature. OCP, N = 2, L = 15. Open symbols, one trial 
move per update. Filled symbols, d = 15 trial moves per 
update. 



VII. OPTIMIZATION 

In Sees. IVl and I VII we presented several ways of updat- 
ing the electric field during charge motion. We also mea- 
sured the mobility of charges. However, the rates Ao have 
been computed in simulation units; time is expressed in 
Monte Carlo trials. As a function of their complexities, 
the different kinds of update require different computa- 
tional effort. In order to choose the best parameters for a 
simulation we should express the efficiency of the various 
versions of the algorithm in terms of CPU time. 

We simulated N = 2 mobile charges in a box of 
size L = 15. T — 0.01 when charges are pointlike 
(temporary spreading) and T = 0.01/ w when spread. 
This set of parameters is representative for simulating 
a monovalent ion in water. At each elementary Monte 
Carlo step (MCS), we try a particle move with prob- 
ability pi = 50%, and a plaquette update with prob- 
ability P2 = 50%. We define a "volume sweep" (VS) 
1 VS = L 3 MCS. 60 000 VS > 2 x 10 8 MCS are per- 
formed after equilibration. Temporary spreading was im- 
plemented in both ways of Sec. I VII first, with steps 1 to 3 
of Fig. ITT1 summed up and stored in a single lookup table; 
second, with multiple steps 2 between each spreading- 
and-recondensing pair of events. 

In Tabled we compare the efficiency of the various up- 
dates introduced in this paper. We used a Pentium 4 at 
2.6 GHz; our C++ code was compiled with an Intel com- 
piler. We conclude that the most efficient field update is 
the temporary spreading of charges on w = 5 cubes. At 
T = 0.01, the mobility reached with w = 5 is close to 
the maximum possible value: D w 0.15 a 2 (sweep) -1 is 
rather close to saturation. We thus do not expect bene- 
fit from further spreading of charges (w ^ 6). As noted 
previously, both versions of temporary spreading yield al- 
most the same mobility. The difference between the two 
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Permanent spreading 


Temporary spreading, 
precomputed current 


Temporary spreading, 
long-ranged particle moves 


A [(VS)- 1 ] 


tcpu [s] 


efficiency [s x ] 


A [(VS)- 1 ] 


tcpu [s] 


efficiency [s 1 ] 


A [VS" 1 ] 


tcpu [s] 


efficiency [s 1 ] 


w — 1 





71 








71 








71 





w = 3 


1.3 x 10" 3 


160 


0.49 


5 x 10~ 2 


592 


5 


6 x 10" 2 


233 


15 


w = 5 


1.4 x 10" 2 


582 


1.4 


1.3 


2372 


33 


1.2 


912 


79 



TABLE I: Comparison of the various algorithms presented. Ao, rate at which field configurations decorrelate, in simulation units. 
£cpu, duration of the 60 000 VS simulation. Efficiency, real rate at which configurations decorrelate, given by Ao x 60 000/£cpu- 
w = 1 stands for the single-link field update, displayed for execution time comparisons. 



is CPU time: long-ranged particle moves lead to a faster 
algorithm thanks to fewer spreading and recondensing 
steps. This version should thus be used for free charges. 

Finally, we have checked that our results remain valid 
for higher densities. We applied our optimal solution 
(temporary spreading over 5 3 sites) to simulate OCPs 
containing N = 14, 34, and 336 positive charges, which, 
respectively, corresponds to number densities n « 0.4%, 
1%, and 10%. 

In Sec. IIIII we calculated a relationship according to 
which Ao oc n\i. This was for Ao in physically relevant 
units of time, like particle sweeps (PS): the effects of each 
charge add up, hence the factor of n. Here we measure 
time in volume sweeps (VS), and work at constant numer- 
ical effort, split amongst particles: the more charges there 
are, the fewer trials each one does. 1 VS = 0.5 /n PS, so 
that A [in (VS) -1 ] = 0.5A [in (PS) _1 ]/ n is directly pro- 
portional to \i. Plotting mobility against temperature in 
Fig-El we find that lowest mobility is found at the lowest 
density; at high density /i decreases, possibly because of 
steric hindrance, but remains greater than when N = 2. 
Thus using our algorithm is always at least as efficient as 
displayed in Table [I] 



VIII. CONCLUSION 



The original version of the local Monte Carlo algorithm 
suffers from two problems at low temperature. First the 
acceptance rate becomes low due to a energy barrier for 
particle motion. A more serious problem is that the mo- 
bility falls even faster than the acceptance rate. We un- 
derstand this fall in mobility by considering the tension 
of the strings left behind particles as they move. The 
different scaling of the acceptance rate and the mobility 
with the spreading parameter w is a clear demonstration 
that two different mechanisms are important in limiting 
particle motion. 

Simple modifications to the algorithm reduce the en- 
ergy barrier for single particle moves, but also the string 
tension. The algorithm is then suitable for simulation 
of lattice models of Coulomb interacting particles. Ex- 
amples include the restricted primitive model for elec- 
trolytes, or lattice models of polyelectrolytes. 




FIG. 14: Mobility of charges temporarily spread on w = 5 
cubes, for an OCP at various densities (n = N/L 3 with L = 
15). 



Combination of the update methods used in this arti- 
cle with the worm update for the transverse field due to 
Alet and S0rensen [l7| can lead to efficient codes for the 
simulation of charge systems at high dilution: Consider a 
set of N charges in a simulation box of size L. It takes a 
computational effort of order NL 2 for these particles to 
diffuse the system size. We have already shown Q that 
the 2L 3 tranverse degrees of freedom of the lattice can 
be integrated over in Oil) sweeps of the worm algorithm 
with an effort scaling as L 3 . One can thus equilibrate a 
dilute system of mobile charges with a computer effort 
which scales as (NL 2 + L 3 ). 



The time moving the particles dominates the time 
needed for the electrostatic integration if N > L, or if 
the density n > 1/L 2 ; when L is large the algorithm re- 
mains efficient even for very dilute charges. It is thus 
well suited to the study of heterogenous systems such as 
surfaces and polyelectrolytes. 
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